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Abstract 

We study some aspects of the invariant pair problem for matrix polynomials, 
as introduced by Betcke and Kressner [3j and by Beyn and Thiimmler |6j. 
Invariant pairs extend the notion of eigenvalue-eigenvector pairs, providing 
a counterpart of invariant subspaces for the nonlinear case. We compute 
formulations for the condition numbers and the backward error for invariant 
pairs and solvents. We then adapt the Sakurai-Sugiura moment method JIj 
to the computation of invariant pairs, including some classes of problems that 
have multiple eigenvalues. Numerical refinement via a variant of Newton’s 
method is also studied. Furthermore, we investigate the relation between the 
matrix solvent problem and the triangularization of matrix polynomials. 

Keywords: matrix polynomials, eigenvalues, invariant pairs, contour 
integral, moments, solvents, triangularization. 


1. Introduction 

Invariant pairs, introduced and analyzed in H, a and [34] . are a gener¬ 
alization of eigenpairs for matrix polynomials. Let P(A) = ^A =0 AA t> e an 
nX7i matrix polynomial, and choose a positive integer k. Then the matrices 
A", S of sizes n x k and k x k, respectively, form an invariant pair of size k 
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for P(X) if 


l 

P(X, S):=J2 A i XSj = °- 

3 =0 

Invariant pairs offer a unified theoretical perspective on the problem of com¬ 
puting several eigenvalue-eigenvector pairs for a given matrix polynomial. 
From a numerical point of view, moreover, the computation of an invariant 
pair tends to be more stable than the computation of single eigenpairs, par¬ 
ticularly in the case of multiple or tightly clustered eigenvalues. Observe that 
the notion of invariant pairs can also be applied to more general nonlinear 
problems, although here we will limit our presentation to matrix polynomials. 

How to compute invariant pairs? Beyn and Thiimmler (0) adopt a 
continuation method of predictor-corrector type. Betcke and Kressner (0), 
on the other hand, establish a correspondence between invariant pairs of a 
given polynomial and of its linearizations. Invariant pairs for P(X) are ex¬ 
tracted from invariant pairs of a linearized form and then refined via Newton’s 
method. 

The approach we take in this paper to compute invariant pairs is based on 
contour integrals. Being able to specify the contour T allows us to select in¬ 
variant pairs that have eigenvalues in a prescribed part of the complex plane. 
Contour integrals play an important role in the definition and computation 
of moments, which form a Hankel matrix pencil yielding the eigenvalues of 
the given matrix polynomial that belong to the prescribed contour. The 
use of Hankel pencils of moment matrices is widespread in several applica¬ 
tions such as control theory, signal processing or shape reconstruction, but 
nonlinear eigenvalue-eigenvector problems can also be tackled through this 
approach, as suggested for instance in [T] and ||4]. E. Polizzi’s FEAST algo¬ 
rithm [31j is also an interesting example of contour-integral based eigensolver 
applied to large-scale electronic structure computations. 

Here we adapt such methods to the computation of invariant pairs. This 
work studies in particular the scalar moment method and its relation with 
the multiplicity structure of the eigenvalues, but it also explores the behavior 
of the block version. We seek to compute invariant pairs while avoiding the 
linearization of P(A), which explains the choice of moment methods. We are 
motivated, in part, by the problem of better understanding the invariant pair 
problem in a symbolic or symbolic-numeric framework, that is, computing 
invariant pairs exactly, or with arbitrary accuracy via an effective combina¬ 
tion of symbolic and numerical techniques: this is one of the reasons why we 
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are interested in the issue of eigenvalue multiplicity. 

The last part of the paper shows an application of our results on invariant 
pairs to the particular case of matrix solvents, that is, to the matrix equation 


P(S) :=i>iS* = 0. 


3=0 

The matrix solvent problem has received remarkable attention in the lit¬ 
erature, since Sylvester’s work [33] in the 1880s. The relation between the 
Riccati and the quadratic matrix equation is highlighted in [7], whereas a 
study on the existence of solvents can be found in [13]. Several works ad¬ 
dress the problem of computing a numerical approximation for the solution of 
the quadratic matrix equation: an approach to compute, when possible, the 
dominant solvent is proposed in S3- Newton’s method and some variations 
are also used to approximate solvents numerically: see for example m, 1221, 
|21j . [27] , The work in [19] uses interval arithmetic to compute an interval 
matrix containing the exact solution to the quadratic matrix equation. For 
the case of the general matrix solvent problem, we can also cite H, h and 
[23]. On the other hand, the question of designing symbolic algorithms for 
computing solvents remains relatively unexplored. Attempts have been made 
to formulate the problem as a system of polynomial equations which can be 
solved via standard methods. However, this approach becomes cumbersome 
for problems of large size (see [2U]). 

Here we exhibit computable formulations for the condition number and 
backward error of the general matrix solvent problem, generalizing existing 
works on the quadratic matrix equation. Moreover, we propose an adaptation 
of the moment method to the computation of solvents. Finally, we build 
on existing work on triangularization of matrix polynomials (see [33] and 
[35]) and explore the relationship between solvents of matrix polynomials in 
general and in triangularized form. 

The paper is organized as follows. Section [2] introduces preliminary no¬ 
tions, definitions and notation. The backward error and condition number 


for the invariant pair problem are computed in Section 2.1 


Section [3] is devoted to the computation of eigenvalues and invariant pairs 
through moments and Hankel pencils. Our main results here consist in The¬ 
orem [5] Corollary [l] and Theorem [7] A comparison of different techniques 


for numerical refinement of invariant pairs is presented in Section 3.4 
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Finally, Sections S@ and [6] are devoted to the applications to matrix 
solvents. 

The methods presented in the paper have been implemented both in a 
symbolic (exact) and in a numerical version. The Maple and Matlab codes 
are available online at the URL 

http://www.unilim.fr/pages_perso/esteban.segura/software.html. 


2. Matrix polynomials and invariant pairs 

A complex n x n matrix polynomial P(A) of degree l takes the form: 

P{ A) = A 0 + AiA + A 2 A 2 + • • • + (1) 

where An ^ 0 and Ai G C nxn , for i — 0,..., i. 

In this work, we assume that P(A) is regular, i.e., detP(A) does not vanish 
identically. 

A crucial property of matrix polynomials is the existence of the Smith 
form (see, e.g., MY- 

Theorem 1. [Thm. Sl.l, uw Every nxn regular matrix polynomial P( A) 
admits the representation 


D( A) = P(A)P(A)P(A), (2) 

where D{ A) = diag (di(A),..., d n (X)) is a diagonal polynomial matrix with 
monic scalar polynomials di{ A) such that di(X) is divisible by d*_i(A); P(A) 
and F(A) are matrix polynomials of size n x n with constant nonzero deter¬ 
minants. 

The polynomial eigenvalue problem (PEP) consists in determining right 
eigenvalue-eigenvector pairs (X,x) G C x C n , with s/0, such that 

P(X)x = 0, 

or left eigenvalue-eigenvector pairs (A, y) G C x C n , with y ^ 0, such that 

y*P( A) = o. 

A particular case of special interest is the quadratic eigenvalue problem (QEP), 
where i = 2. Typical applications of the QEP include the vibrational anal¬ 
ysis of various physical systems. A considerable amount of work has been 
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done on the theoretical and computational study of the QEP: see for instance 

Id- 

As for the linear case, there is a notion of algebraic and geometric multi¬ 
plicity of eigenvalues of matrix polynomials. If A 0 is an eigenvalue of P(A), 
the algebraic multiplicity of A 0 is its multiplicity as root of detP(A), whereas 
the geometric multiplicity of Ao is the dimension of the null space of the 
matrix P(Ao). 

Invariant pairs , introduced and analyzed in |3j and |6] , are a generaliza¬ 
tion of the notion of eigenpair for matrix polynomials. 

Definition 1. A pair (A, S') G C nxk x C kxk , X ^ 0, is called an invariant 
pair for P(A) if it satisfies the relation: 

P(X, S ) := A e XS e + • • • + A 2 XS 2 + AiXS + A 0 X = 0, (3) 


where Ai G C nxn , i = 0,and k is an integer between 1 and In. 


The following definitions proposed in |3| and D3| will be helpful for our 
work, for instance, to allow for rank deficiencies in A". 


Definition 2. A pair (A", S ) 
m G N such that: 


G C nxk x C kxk is called minimal if there is 
~X S'™” 1 ' 


V m {X,S) 


X s 
X 


has full column rank. The smallest such m is called minimality index of 

(X,S). 


Definition 3. An invariant pair (A", S) for a regular matrix polynomial P(A) 
of degree I is called simple if (A, S) is minimal and the algebraic multiplicities 
of the eigenvalues of S are identical to the algebraic multiplicities of the 
corresponding eigenvalues of P( A). 


Invariant pairs are closely related to the theory of standard pairs pre¬ 
sented in HZ], and in particular to Jordan pairs. If (A", S) is a simple invariant 
pair and S is in Jordan form, then (A, S) is a Jordan pair. 
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As an example, consider the following quadratic matrix polynomial, dis¬ 
cussed in B3: 


SA 

to 

'1 0 o' 
0 1 0 

+ A 

"-2 0 1' 
0 0 0 

+ 

'1 0 0" 
0-10 


0 0 0 


0 0 0 


0 0 1 


(4) 


It has eigenvalues Ai = 1 with algebraic multiplicity 3 and A 2 = — 1 with 
algebraic multiplicity 1. A corresponding Jordan pair (.X, J ) is given by: 


X 


"0 

0 

1 

o' 

/ 

1 — 

1 

1 

0 

1 

, J = diag -1,1, 

1 1 

0 1 

0 

0 

0 

0 

V 


The notion of invariant pairs offers a theoretical perspective and a numeri¬ 
cally more stable approach to the task of computing several eigenpairs of a 
matrix polynomial. Indeed, this problem is typically ill-conditioned in pres¬ 
ence of multiple or nearly multiple eigenvalues, whereas the corresponding 
invariant pair formulation may have better stability properties. 

I 11 particular, simple invariant pairs play an important role when using 
a linearization approach as in [3], and ensure local quadratic convergence of 
Newton’s method, as shown in [|23|; see also [53] . 


2.1. Condition number and backward error of the invariant pair problem 
In the following sections, we present explicit formulas for the backward 
error and the condition number of an invariant pair (A", S'). We follow the 
ideas presented in the articles [36] and [21] , which give expressions for back¬ 
ward errors and condition numbers for the polynomial eigenvalue problem 
and for a solvent of the quadratic matrix equation. 


2.1.1. Condition number 

Let (A, S') be an invariant pair for P(A), and consider the perturbed 
polynomial 

(P + AP)(A) = (Aq + AA 0 ) + A(Ai + AAi) + • • • + A + A Af), 


1 

where AP(A) = X l AA i and AA 0 ,..., A £ C nxn . Let (AX, A S') be a 

i=0 

perturbation of (X, S') such that 


(P + AP)(X + AX, S + AS) = 0. 
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A normwise condition number of the invariant pair (A", S) can be defined as: 


k(X, S ) = limsup < 


AA 

AS 


X 

s 


A : (P + AP)(X + AX, S + AS) = 0, 

||AAj||ir < eai,i = 0,... (5) 


The cti are nonnegative weights that provide flexibility in how the pertur¬ 
bations are measured. A common choice is cc* = ||A,||f; however, if some 
coefficients are to be left unperturbed, A A* can be forced to zero by setting 

ai = 0 . 


Theorem 2. The normwise condition number of the simple invariant pair 
(A, S ) is given by: 


where 


k{X, S ) 


[B x -Bs] + Ba 



( 6 ) 


t 


b *=E([( 5 T®A]), 

j=o 

B a = [ai(XS e ) T <g) I n •• 


= ((^■ l_1 ) T ® AjXST ), 

3 =1 i= 0 

a 0 X T <g> I n ] . 


Proof. By expanding the hrst constraint in (|5]) and keeping only the first 
order terms, we get: 


t i l 

AAjXS j + AjAXS j + Yl AjXBS j (AS) = 0(e 2 ), 

3=0 3=0 3=1 


( 7 ) 


where D S 3 denotes the Frechet derivative of the map S 3 S 3 : 

j -1 

: AS ^^FASS*-*- 1 . 

i =0 
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Using on equation 0 the vec operator and its property (see [IS] , pp. 28): 

vec (AXB) = ( B t ® A) vec (X), (8) 


we obtain: 


»vec(A P(X, S)) = vec E AA > xsj = E vec = E ([( X S :l ) r ® J„1 vec (A/1 ;/ j) = 

2=0 / 2=0 2=0 
'vec (A Af)/af 


= [a l {XS l ) T ®I n ••• a 0 X T (g>/„] 


vec (AA 0 )/ao- 


=: Ba vec (AA), 


(l \ l l 

E AjAXS j = Evec(A,AA'Sb = E ([(5^) T ® A,-]) vec (AX) 
2=0 / 2=0 2=0 


=: -Bx vec (AX), 


2=1 

2-1 


2-1 


vec j2 A j xosj ( As ) = vec J2 A i x Yl siAssj - i ~ 1 =EE vec ( i i IS ' ASSi "" 1 ) = 


2=1 *=0 


l j-1 


2=1 i= 0 


= E E ((‘S ,j_i_1 ) T ® AjXS*) vec (AS) =: vec (AS). 

[B x B s \ y = —B a x + 0(e 2 ), 


2 = 1 1=0 

Then, we have: 


where 





vec (A A e )/a e 

y = 

vec (AX)' 
vec (AS 1 ) 

, and x = 




vec (A A 0 )/a 0 


and therefore 


Wvh 


Tvec (AX)' 


" AX 1 

vec (AS 1 ) 

2 

AS \ 


So we have that the definition (J5]) is equivalent to the following 


lim sup 
£->0 


1 lll/lla 


e 

' X ' 



_ s _ 

F 


[-Bx B s ] y 


-B A x + 0(e 2 ),\\x\\ 2 < e 


[B x 

B s ] + B A 

2 


' X ' 
S _ 

F 


where the matrix [Bx B$] has full rank if the invariant pair (A", S) is 
simple (see [Thm. 7, 0]). □ 






































In order to better illustrate Theorem[2j let us consider the particular case 
k — 1. When k — 1, invariant pairs (XfS) coincide with eigenpairs (a;, A). 
In this case, the matrices Bx, B$ and B A in (|6| 


are: 


B x =E (t( AJ ) T ® A A) = E ( Ai ^) = p ( A )- 

j'=o j=o 

ss=EE(( Ai " < " 1 ) T ®^* A< ) = EE( aH v) =^( a )*. 

j— 1i—0 j=li=0 

T ^ 


= [aeX e x T 0 I n ac-i\ e 1 x T 0 /, 

Note that: 


aox 


Bax = \aiX £ x T > 


a 0 x T (g) I n 1 


vec (AA e )/a e 


vec (Aj4o)/ao 


= vec + • • • + ATox) = vec (AP(X)x) = AP(X)x. 

Therefore, we obtain: 


[Bx B S \ y = ~B a x + 0(e 2 ) 4^ [P(A) P'(A)x] 
4^ P(A)Ax + P'(X)xAX + AP(A)x = 0(e 2 ). 


Ax 

AA 


= —AP(A)x + 0(e 2 ) 


The last equation is consistent with the first part of the computation of 
the condition number for a nonzero simple eigenvalue A of P(A) presented 
in [Thm. 5, [35]]. The second part differs, because here we are estimating 

/\ T, 

, whereas classical condition numbers for eigenvalue problems 

J f 

typically take into account angles between left and right eigenvectors. Of 
course, it would also be interesting to formalize a similar approach for in¬ 
variant pairs, based on angles between suitable matrix manifolds (such as 
partially developed in [3]). 

2.1.2. Backward error for P(X, S ) 


Let ttj, for i = 0, ...,£, be nonnegative weights as in Section 2.1.1 The 
backward error of a computed solution (. X , S) € C nxfc x <C kxk to (J3]) can be 
defined as: 


V (X,S) = min{e: (P + AP)(X,S) = 0, || AA^p < ea u i = 0,... ,£} (9) 
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( 10 ) 


By expanding the first constraint in (J9]) we get: 

- P(X, S) = A A £ XS e + • • • + AA 0 X. 

Then, we have 

-P(X,S) = [aJ 1 AA e ... a^AA! a^AHo] 


'a t XS t 

ai XS 
CtQ X 


Taking the Frobenius norm, we obtain the lower bound for the backward 
error: 

„(.Y ,§) > __„_. 

(o? t \\XS% + ■ ■ ■ + a?||A'S||J, + ag|| X|| J.) 1 / 2 

Consider again equation © Using ((8|) , we obtain: 

- vec (P(X, S)) = ((XSY 0 I n ) vec (A At) + • • • + (X T 0 I n ) vec (AA) 

'vec (A A t )/ai 

= [a e (XS e ) T 0l n ... a 0 X T 0 I n ] 


vec (AA 0 )/a 0m 


which can be written as: 


Hz = r, He C nfcx(/+1 >* 2 (11) 

Here we assume that H is full rank , to guarantee that © has a solution 
(backward error is finite). Then the backward error is the minimum 2-norm 
solution to: 

0 (JY,S) = ||H+r|| 2 , (12) 

where H + denotes the Moore-Penrose pseudoinverse of H + . 

Eq. (12) yields an upper bound for rj(X, S): 


r](X,S) < ||tf + || 2 ||r|| 2 = 


\T 2 


^min (H) 


where cr m i n denotes the smallest singular value, which is nonzero by assump¬ 
tion. Note that: 


<7min (H ) 2 = A mia (HH*) = A min (aJ(XS e ) T XS e »/„ + ••• + a\X T X ® I n ) > 
> ag(j m i n (X + • • • + a^cr m i n (AS ') 2 + aQ<T m i n (A) 2 . 
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Thus we obtain the upper bound for r](X, S): 


§)< _ ll^,g)llF _ 

’ ~ (aja^XS*)* + • • • + aja^XSy + a 2 0 a min (X) 2 

In the particular case k = 1, the approximate invariant pair (. X , S ) co¬ 
incides with an approximate eigenpair (x, A). In this case, the definition (J9]) 
becomes: 



ri(X 


r)(x, A) = min{e : (P + A P)(x, A) = 0, || AAj||p < ect;, i = 0,..., £}, 

which is the definition of the normwise backward error of an approximate 
eigenpair (x, A) for P( X)x = 0, presented in [(2.2), plb] ]. 

2.1.3. A numerical example 

Let us consider the power plant problem presented in [5] and in |32]. This 
is a real symmetric QEP, with P(A) of size 8 x 8 , which describes the dynamic 
behaviour of a nuclear power plant simplified into an eight-degrees-of-freedom 
system. The problem is ill-conditioned due to the bad scaling of the matrix 
coefficients. 

The maximum condition number for the eigenvalues of P(A), computed by 
the MATLAB function polyeig, is: 


Kmax = max condeig A = 1.0086e+008. 
AeA 


Using the method that will be presented in Section |3.2| and Section |3.4.1 
we compute an invariant pair (X, S) associated with the 11 eigenvalues with 
largest condition number inside the contour T = 7 + pe lB (7 = 80 + 10 *, 
p = 170). The condition number and backward error for (.X , S) are 

k(X, S ) = 565.6746 and p{X , S) = 4.4548e - 017. 

Observe that k(X, S ) is significantly smaller than K max . 


3. Computation of invariant pairs 

Numerical methods based on contour integrals for the computation of 
eigenvalues of matrix polynomials and analytic matrix-valued functions have 
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recently met with growing interest. Such techniques are related to the well- 
known method of moments, where the moments are computed by numerical 
quadrature. 

In this section we explore a similar approach for computing invariant 
pairs. Our main reference is the Sakurai-Sugiura method (see [T] and |32j). 
as well as the presentation given in [3]. 

3.1. The moment method and eigenvalues 

Let us begin by briefly recalling a few basic facts about the Sakurai- 
Sugiura moment method. Here we essentially follow the presentation given 
in Hj. 

Let T be a closed contour in the complex plane and let u and v be arbi¬ 
trarily given vectors in C n . Define the function: 

/(A) := u H P( A)-V 

In the following, it will be understood that no eigenvalue of P( A) should lie 
exactly on the contour T : each eigenvalue should be either inside or outside 
the contour. 

The next theorem, which can be found in |[T], gives a representation for 
/(A) that will prove useful later on. 

Theorem 3. [Thm. 3.1, Jl]] Let D( A) = diag (rfi(A),..., rf n (A)) be the 
Smith form of P( A), and let E( A) and F( A) be as in (|2l). Let yJX) = 
u H q j (X)p j (X) H v, 1 <j< n. Then 


(13) 

j =i J 

where q,(A) and p-(A) are the column vectors of E( A) and F(X) H , respec¬ 
tively. 

Definition 4. Let fceN. The k-th moment of f(z) is: 

^ k= 2ln?f zk -f^ dz - 
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For a positive integer m, define the Hankel matrices H 0 , Hi e C mxm as 
follows: 


ho 

hi 

l^m— 1 


hi 

h2 

t^m 

hi 

h2 ' 


, Hi = 

h2 

h3 •• 

/^m+1 

hro-1 

hm 

‘ /^2m—2 


Pm 

hm+1 

l^2m— 1 


The eigenvalue algorithm presented in [1] relies on the following result: 

Theorem 4. [Thm. 3.4, Suppose that the polynomial P( A) has exactly 
m eigenvalues Ai,..., A m in the interior ofY, and that these eigenvalues are 
distinct, simple and non degenerated. If Xn{Xf) 7^ 0 for 1 < I < m, then the 
eigenvalues of the pencil Hi — A H 0 are given by Ai,..., A m . 

We now wish to investigate the behavior of the above method when the 
hypothesis that the Aj’s are distinct is removed. In particular, we aim to 
generalize Theorem [4j which is based on the Vandermonde factorization of 
//(i and Hi. 


Theorem 5. Suppose that P(A) has exactly m eigenvalues in the interior 
ofT, namely, distinct eigenvalues A 0 ,..., X s with algebraic multiplicities m 0 , 

..., m s , respectively, such that m = m 0 H- Vm s . Moreover, assume that no 

eigenvalue of P{ A) lies exactly on the contour T. If the geometric multiplicity 
of the A i’s, for i = 0,..., s, is equal to one, then the matrix H 0 is nonsingular 
and eigenvalues of the pencil Hi — A H 0 are given by X 0 ,..., X s with algebraic 
multiplicities m 0 ,..., m s . 

Proof. 

Suppose that, in the Smith form Q of P(A), the matrix D( A) is in the 
form D{ A) = diag(rfi(A),..., d n _i(A), d n ( A)), where 

r 

4(A) = (A - A 0 ) mo (A - AO " 11 ■ • • (A - A s ) m ‘ (X - A,) m % 

i=s+l 

and A s +i,...,A r are the eigenvalues of P(A) located outside the contour 
T, with algebraic multiplicities m s+ i,... ,jn r . Moreover, define d n ( A) = 
n: =0 (A — Aj) m % i.e., d n (X) is the factor of d n (X) whose roots are the eigen¬ 
values of P(A) located inside T. 
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Since the geometric multiplicities of the A,’s are all equal to one, the 
factors (A—A*), for i — 0,..., s, do not appear in the monic scalar polynomials 
^o(A), • • * 5 d n — l(A). 

Applying Theorem [3j we have: 



where 

n j / \ 

<p{z) ='52xj(z)h j (z), with hj(z) = 

j = i a j\ z ) 

We can introduce partial fraction decompositions and write 



where c ]:i G C, for j = 0,..., s and i = 1 ,...,rrij. Classical results on 
residues then yield 


S ™j 

dk EE c^jRes 

3=0 i= 1 


0 - A j) 


^3 = 


g TILj ■ ^ 

EE c »7iWT 1 - A 

3=0 i=l 

s i . dr 1 


(i — 1)! \j dz l ~ l \ (z — A j)' 1 


EE^y^iyr^z^iC) 

j=0 i=l x 7 
s m j 


i k—i+1 

S 


(16) 


3=0 i =1 
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where 


^j,i 


-^y (k - i + 2)0 - i + 3) • • • k 
0 otherwise. 


if 


k > i — 1 , 


Now, consider the pencil Hi—XH 0 with H 0 and H i defined as in (15). Because 
of (16), and of the fact that Aq, .. ■, A s are roots of d n ( A), the moments /i^ 


satisfy a linear recurrence equation of the form: 


k — dm-l^k-l + + ' ' ' + Oo/ifc —m* (17) 

Moreover, d n ( A) is the polynomial of smallest degree that has roots Ao,..., X s 
with the prescribed multiplicities mo, • • • , m s , so the recurrence 0 has the 
shortest possible length; also, note that the ads in © are actually the 
coefficients of d n ( A). Therefore, the matrices 77 0 and H\ have full rank. The 
same argument shows that H 0 and H\ are rank-deficient if taken of size larger 
than m x m (see |29], [[T 6 ], Vol. 2 , pp. 205)] and | 8 j). 

As a consequence of the shifted Hankel form of Ho and H i, we have 
7/(i 6' = 77i, where C is a matrix in companion form 


"o 

0 

... 0 

x 0 

1 

0 

... 0 

Xi 

0 

1 

... 0 


0 

0 

... 1 

%m—l 


and its last column is given by the solution of the linear system 



Xo 


/^m 

H 0 

Xi 

= 

Mm+1 


%m— 1 


l^2m—l 


The polynomial of degree m: 

p(\) =\ m - - X 0 


(18) 


is a scalar multiple of (A — A 0 ) mo (A — Ai) mi • • • (A — A s ) ms , and its roots are 
the Aj’s generating the entries of the pencil H\ — A77 0 . So we also have that 
the /ids satisfy the recurrence 0 : 


H'k ■^m—l/^k—1 ^m—2f^k—2 * * * “h Z'OH'k— m? 
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where k = m, m + 1 ,..., 2m. 

Consider now the Jordan matrix 

'Ji 
J = 

J s_ 

where each block Jj, of dimension rn l . is a square matrix of the form 

'A,: 1 

T _ Aj 

Ji i 

l 

A,. 

and define the confluent Vandermonde matrix 

V = (v J T v • ■ ■ (J T ) r " 1 v) , 

where v 7 = ^e% ni ^ T • • • j s partitioned conformally with J and 

e'”' ,r =(l 0 ■■■ 0 ) T is the rrq-dimensional unit coordinate vector. Then 

we have 

VC = (v J T v ■ • • ( J T ) r ~ 1 v) C 

= (J T v (J T ) r-1 v -(x 0 / + xiJH -f x m _i J r ~ l ) T v) 

= (J T v ■ ■ ■ ( J T ) r ~ 1 v (J T ) r v) 

= J T v, 


where we have used the property p(J) = 0, that is, the Cayley-Hamilton 
theorem. 

We can now introduce the Vandermonde decomposition of the Hankel 
matrices Hq and H\. From the results presented in [ 8 ] it follows that there 
exist block matrices B 0 = diag(-Dj°\ ..., D i°^) and B\ = diag(Ij| 1 \ ..., D^), 
partitioned conformally with J, satisfying the conditions BqJ 1 = JBq and 
B\J 7 = JBi, so that 


Hi = V T BiV , for i = 0 , 1 . 
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Moreover, we can prove that JB 0 = Bp. 

HqC = H 1 ^ ( V t B 0 V)C = V t B 1 V V t B 0 J t V = V t B 1 V 
V t JB 0 V = V t B 1 V JB 0 = B u 

where we used the properties VC = J T V and B, L J 7 = JBi. 
Therefore, we have: 

H x - XH 0 = V t B 1 V - A V t B 0 V = V t JB 0 V - XV T B 0 V 
= V T (J - XI)BqV. 


, A s with respective multiplicities 

□ 


So the eigenvalues of H i — A H 0 are Ao, 
m 0 ,.. .,m s . 

Example 1 . Consider the matrix polynomial: 


P{ A) = A 2 / + A 


P( A) has eigenvalues: X\ = ^,A 2 = 1 ,A 3 = 2, A 4 = 3 with algebraic multi¬ 
plicities: m 1 = 2 , m2 = 2, m3 = l,m4 = 3. The associated Smith form is: 


that: 


"-2 

0 

0 

0 ' 


'1 

0 

0 

1 ' 

0 

-1 

0 

0 


0 

1 

0 

0 

0 

0 

-6 

0 

+ 

0 

4 

0 

9 

0 

0 

0 

0 

-5 


0 

0 

0 

6 


Hn = 


Then, we have: 


-3 

-7 

-9 

-7 

-9 

-21 

2 

-9 

-21 

2 

-12 

-21 

-12 

-109 

2 

8 


!) 2 (A- 

-1) 2 (A-2)(A 

~ 3) 3 ) 



[2 -2 

1 — l] and v = 

' o' 

O 

T 

, we 

-21 - 
2 


'-7 

_q =21 
v 2 

-12' 

-12 

, H,= 

-9 

tF -12 

-109 

-109 

-21 

12 -f 9 

8 

-123 

8 

-123 

8 J 


2 

-12 

-109 -123 

8 8 

8 

-551 
32 J 


C = 


0 0 0 ^' 

1 0 0 | 

1 0 


0 

0 0 1 


Note that the eigenvalues of C are 1 , 1 . Moreover, the companion matrix 
C is associated with the monic polynomial: 

P( A) = A 4 


3A 3 + -A 2 
4 


3 1 

2 A+ 4 ’ 


17 



whose roots are indeed 1,1. 

/n fact, the Vandermonde factorization for H 0 and Hi is Hi = V T BiV, 
i = 0,1, where 


'1 

1 

2 

1 

4 

1 “ 



' 0 

-2 

0 

0 ' 


-2 

-1 

0 

0 ' 

0 

1 

1 

! 


B 0 = 

-2 

0 

0 

0 

, Bi = 

-1 

0 

0 

0 

1 

1 

1 

4 

1 


0 

0 

—3 

—2 

0 

0 

—5 

—2 

0 

1 

2 

3 



0 

0 

-2 

0 


0 

0 

—2 

0 


and JB 0 = £> 1; with 

1 1 0 O' 

t 0^00 

J= 0 0 1 1 ' 

0 0 0 1 

Example 2. Consider the quadratic matrix polynomial: 


P( A) = A 2 

' 1 0 O' 

2 10 

+ A 

" 0 0 o' 

-4 -2 0 

+ 

"-2 1 -2' 
2 10 


1 

1 

to 


2-2 4 


-1 1 -2 


with associated Smith form: 


D(A) = diag((d 1 (A),d 2 (A),d 3 (A)) 


1 0 0 

0 (A — l) 2 0 

0 0 (A — 1) 3 (A + 1) 


The Jordan matrix associated with the linearized form of P( A) is: 

j _ J\ 0 
[o J 2 \ ’ 

where each block Ji is: 

-1 0 0 O' 

0 110 

0 0 11 ’ 

0 0 0 1 


Jl = 


1 1 
0 1 


J 2 = 


Note that we have the eigenvalue A = 1 in different Jordan blocks. 

Choose T as the circle <p(t) = 1 + which contains 5 eigenvalues A = l, 

and consider vectors u = [ 3 1 — 2] T and v — [3 —1 — 2] T . Theorem\5 
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implies that the moment method yields the eigenvalues of P( A) inside the 
contour, which are roots of d 3 (\), i.e., Ao = 1 with multiplicity mo = 3. Let 
us compute the matrices Ho and Hi of size 3x3: 


"7 

3 

3' 


'3 

3 

7' 

3 

3 

7 

, H 1 = 

3 

7 

15 

3 

7 

15 


7 

15 

27 


The matrix Ho is nonsingular. Then, we have: 


C 


H^Hi = 


0 0 1 
1 0 -3 
0 1 3 


Note that eig(C) = 1,1,1. 

As the contour contains 5 eigenvalues, we might ask what happens when 
taking the Hankel matrix Ho of size 5x5: 


H 0 = 


7 3 3 7 15 

3 3 7 15 27 

3 7 15 27 43 

7 15 27 43 63 

15 27 43 63 63 


This matrix is singular, therefore we will not be able to find all the 5 eigen¬ 
values inside the contour. The method miss the additional multiplicities as¬ 
sociated with the polynomial d 2 ( A). 

Remark 1. We conclude that the scalar moment method can be used to 
compute the (possibly multiple) eigenvalues of P( A) that belong to the interior 
of T and that are roots of the polynomial d n ( A). The method misses the 
additional multiplicities associated with the polynomials di(X),..., d n _ i(A). 


The above remark is consistent with the fact that the Jordan form of a 
companion matrix only contains one Jordan block for each eigenvalue: it is 
not possible to capture multiple eigenvalues associated with several Jordan 
blocks. 

In order to “see” the additional eigenvalues that are roots of di(A),..., d n - 1 
the block version of the moment method may be useful (see Section 3.3). 


(A), 
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3.2. Computing invariant pairs via moment pencils 

Let T be a closed contour, Ai,..., A m all the eigenvalues of P(A) in the 


interior of T and the matrices H 0 and H i defined as in (15). 

For k — 0,1,..., m — 1 and a nonzero vector v e C n , consider the vectors: 


Sk = 2^,% Z P(Z) ~ Viz ' 


(19) 


The method proposed in [1] for the computation of the eigenvectors of P(A) 
is based on the following result. 


Theorem 6 . [Thm. 3.5, mi Let (Ai, Wi), i — 1,..., m be eigenpairs for the 
matrix pencil Hi —XH 0 , where the simple, distinct, nondegenerate eigenvalues 
Aj belong to the interior of a given closed contour T. Let S = [s 0 , ■ ■ ■ , s m -i]- 
Then, for i — 1,..., m, the vector r/j = Swi is an eigenvector of P(A) corre¬ 
sponding to the eigenvalue A*. 

Theorem [6] is readily applied to invariant pairs. 

Corollary 1. With the hypotheses of Theorem [fij, S = [so, 
and C = Hf ] H\. Then the pair (S, C) satisfies P(S,C ) = 
a simple invariant pair for P( A). 

Proof. Note that the pair (Y, A), where A = diag(Ai,..., A m ) and Y = 
[yi,... ,y m ], is clearly an invariant pair for P(A), that is, 


Si, • • ■; s m _ij 
0 , i.e., (S,C) is 


i 

P(Y, A) = A i YAj = °- 

3=0 


Moreover, we know that C = Hf l H\ = V~ l AV , where V is the classical 
Vandermonde matrix associated with Ai,...,A m , and that the columns of 
V -1 are eigenvectors of H\ — A H 0 . So we have 


0 


E A i YA> = E A i YA>v 


3=0 


3=0 


t t 

AjYVV- x h?V = J2 A 3 SC 3 = P{S,C), 


3=0 


3=0 


that is, (S,C) is also an invariant pair of P(A). 


□' 
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What can we say about more general cases, where some of the hypotheses 
of Theorem [ 6 ] are removed? If we remove the hypothesis that the Ads are 
distinct, we can prove the following. 

Theorem 7. With the hypotheses of Theorem^ let S = [s 0 , Si, ..., s m -i] 
and C = H^Hi. Then the pair (S', C) satisfies P(S,C ) = 0, i.e., ( S,C ) is 
a simple invariant pair for P( A). 


Proof. Consider again the columns qi(X ),..., q n ( A) of the matrix F( A) in the 
Smith form (| 2 ]) and the dehnition of sy 0 given in (19). A similar computation 
to (16) shows that 

S = [s 0 , • • •, s m _i] = QV, 

where 


Q — [Qo, ■ ■ ■, Qs[i 

Qj = [7oj?n(Aj), 7i• • • > 7m i -ij?i mj '“ 1) (Aj)], for 3 = 0, ..., s, 

the yjj’s are complex coefficients and V is the confluent Vandermonde matrix 
defined above. 

It is shown in [1] (Lemma 2.4) that, if a complex number £ is a root of 
dj( A) for some index 1 < j < n, then -P(C)< 2 j(C) — 0 . In our case, this implies 
that the vector q n ( A) is a root polynomial of P( A) corresponding to the 
eigenvalue \ 3 . for each j = 0 ,..., s; see ra, section 1.5, for the definition and 
properties of root polynomials. It follows that [q n (Xj), q' n (Xj), ..., qi 1 ^ A, ' ) ] 
forms a Jordan chain for the eigenvalue A j. So we have that (Q, J) is an 
invariant pair for P( A). Moreover, if C — Hq 1 Hi as usual, we have 


0 


AjQj’ = Yi A iQ J ’ v = 


3 =0 


3=0 


l l 

Y^AjQVV-'J’V = J2 A jSC j = P{S,C), 


3=0 


3=0 


Therefore, (S,C) is a simple invariant pair for P(X). 
Example 3. Consider the matrix polynomial: 

1 O' 

0 0 ’ 


P( A) = A 2 


1 0 
0 1 


+ A 


-2 0 


+ 


□ 
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which has eigenvalues \\ = 0 with algebraic multiplicity 1 and A 2 = 1 with 
algebraic multiplicity 3. 

Suppose we are interested in the eigenvalues X 2 . Then we can choose a con¬ 
tour r(t) = z 0 + Re lt , t G [0, 2tt], where z 0 — 1 and R = 

Choosing the vectors u — [l — l] T and v = [—1 l] T ; we find that: 



'-1 -2 

-5' 


' -2 

-5 

- 10 ' 

H 0 = 

-2 -5 

-5 -10 

1 1 

° 

, H x = 

-5 

-10 

-10 

-17 

-17 

-26 


Then, we have that the pair ( S,C) given by Theorem [?| 


S = 

0 -1 - 2 ' 
113 

and C = H^H l = 

"0 0 1 ' 
1 0 -3 




0 13 


is an invariant pair, i.e., it satisfies P(S,C ) = 0 . 

Note that the companion matrix C is associated with the monic polynomial: 

p{ A) = A 3 -3A 2 + 3A-1, 
which has as roots: 1,1,1. 

3.3. The block moment method 

Instead of the scalar version of the moment method, we can consider 
a Hankel pencil constructed by block moments M k € C^ x ^, for a suitable 
positive integer £. 

Definition 5. Let k be a positive integer and U, V G C” x * nonzero matrices 
with linearly independent columns. For k = 0,1 ,..define the block moment 
M k G as: 

M k = ^~ l z k U H P(z)~ 1 Vdz. 

2m Jy 

Then the block Hankel matrices H^ 0 , H^i G C m ^ xm ^ are defined as: 



" M 0 

Mi ■ 

Mfn-l 


'Mi 

m 2 ■ 

• Mfn ' 

^ o = 

M, 

m 2 ■ 

■ Mfn 

, H (1 = 

m 2 

m 3 • 

Mfh+i 


Mfn-l 

Mfa • 

M 2 m-2 


Mfh. 

Mfh+i ■ 

M 2 fh-i 
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Polynomial eigenvalue computation via the eigenvalues of the pencil H^i — 
XH^o is discussed in [lj and [4]. See also [26J for an application to acoustic 
nonlinear eigenvalue problems. 

Invariant pairs can be computed from block moments by applying an ap¬ 
proach that is similar to the one described in the previous section for the 
scalar version. For k = 0,1,..., m — 1, consider the matrices Sk G C nx ^ 
defined as: 

s ‘=2 H ztp ^~ 1Vdz - 

Then, we have the following result. 


Proposition 1. Let V be a closed contour, let the block Hankel matrix 0 € 
C mixnne, non singular an d m the number of eigenvalues inside ofY. If 
rrif = m and Y = [So, • • •, S^-i], T = Hf (j l Hg\, then the pair (Y, T ) satisfies 
P(Y,T) = 0, i.e., (' Y,T) is a simple invariant pair for P( A). 

With the condition that the size of the block Hankel matrix H^ 0 G C m ^ xm ^ 
is equal to the number of eigenvalues inside of P, i.e., if rrif = m, we get: 



"o 

0 

... 0 

1 

1 


/ 

0 

... 0 

-Ad 

T = H-jHs i = 

0 

I 

... 0 

<M 

^ •• 
1 


o • 

• o 

■ ■ I 

Xm— 1 


where 


' -Ad ' 


" M m ' 

-Ad 

= HI o 1 

M m+ i 

i 

i 

£ 

1 


A?2m- 1 


Consequently, since T has a block companion form, the problem of find¬ 
ing the eigenvalues Ai,..., A m of is equivalent to the problem of finding the 
eigenvalues of the matrix polynomial: 


L{ A) := X 1 + AViA £_1 + • • • + XiA + X 0 = 0. 
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Example 4. Consider again the matrix polynomial of Example^ 


P( A) = A 2 


' 1 0 o' 


" 0 0 O' 


"-2 1 -2' 

2 10 

+ A 

-4 -2 0 

+ 

2 10 

-1 1 -2 


2-2 4 


-1 1 -2 


with associated Smith form: 


Z>(A) = diag((d 1 (A),d 2 (A),d 3 (A)) 


1 0 0 

0 (A — l) 2 0 

0 0 (A — 1) 3 (A + 1) 


In Example [|, we found that the scalar moment method, i.e. when f = l, 
missed the additional multiplicities associated with the polynomial d 2 ( A). 
Consider now £ = 2 as the size of the block moments M k , the contour ip(t) = 
1 + jqC 1 , containing 5 eigenvalues A = 1, as before, and the matrices: 



1 0 


1 3 

U = 

5 -3 
2 -4 

, v^ = 

0 1 
-2 4 


We find the block moments M k : 


M 0 = 
M 3 = 


-9 -12 
9 12 

-21 30 

15 -15 


M 1 = 


-1 -22 
-1 27 


Mo = 


-5 -8 
1 18 


, Ma = 

'-49 

92 ' 

, M; = 

'-89 

178 ' 


41 

-72 

1 0 

79 

-153 


Then, we have the Hankel matrix H L0 : 


M 0 

M 1 

m 2 

Mi 

m 2 

m 3 

m 2 

m 3 

m 4 


-9 

-12 

-1 

9 

12 

-1 

-1 

-22 

-5 

-1 

27 

1 

-5 

-8 

-21 

1 

18 

15 


-22 

-5 

-8 

27 

1 

18 

-8 

-21 

30 

18 

15 

-15 

30 

-49 

92 

-15 

41 

-72 


The matrix 0 is singular. This happens because there are just 5 eigenvalues 
inside the contour and H^q has size 6x6. Then, we have to reduce the 
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matrices H^ 0 and H^i to match the number of eigenvalues in the contour. 
Therefore, we get the truncated matrices: 


'-9 

-12 

-1 

-22 

-5' 


'-1 

-22 

-5 

-8 

-21' 

9 

12 

-1 

27 

1 


-1 

27 

1 

18 

15 

-1 

-22 

-5 

-8 

-21 

(Try 

II 

-5 

-8 

-21 

30 

-49 

-1 

27 

1 

18 

15 

1 

18 

15 

-15 

41 

—5 

-8 

-21 

30 

-49 


-21 

30 

-49 

92 

-89 


Then, we obtain: 


T = = 


0 

0 

0 

-2 

1 

0 

0 

0 

-1 

0 

1 

0 

0 

4 

-3 

0 

1 

0 

2 

0 

0 

0 

1 

-2 

3 


The eigenvalues of the matrix T are 1,1,1,1, l, which are all the eigenvalues 
inside the contour. 

Moreover, computing the matrix Y = [5*0, Si, 5*2], using we get: 


Y 


0 1 
0 -2 



2 0 
0 0 
-3 -4 


Then, ( Y,T ) is an invariant pair for P(X). 

Experimentally, we noted that the block method allows us to better “cap¬ 
ture” the multiplicity structure of eigenvalues, when there are several Jordan 
blocks per eigenvalue. Further investigation of this approach will be the topic 
of future work. It should be pointed out that the results in [3], and partic¬ 
ularly Theorem 3.3, provide useful insight into a generalized block moment 
method and into the (good) behavior of the method in presence of multiple 
eigenvalues. 

Another delicate issue pertaining to contour integral method is the choice 
of T. If some information about the localization of the eigenvalues is available, 
one can choose the contour accordingly. In other cases, T may be taken as a 
circle for ease of computation, as we do here. 

A related question is: how many eigenvalues of P(X) live inside a given 
contour? Even an approximate estimate can be useful to choose T and k 
consistently. An answer to this problem is provided in HD and [15] . In 
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particular, Theorem 2 in [IT] points out that the number m of eigenvalues of 
P(A) that are inside T is given by: 


m = 


y trace 



( 20 ) 


For practical computation, the right hand side of equation (20) can be ap¬ 
proximated by a quadrature rule, thus yielding an estimate for m. 

Moreover, the choice of T can be combined with shifting techniques for 
the eigenvalues of P(A): see for instance 


3-4- Numerical approximation and refinement of invariant pairs 

When implementing numerical computation of invariant pairs via the 
(block) moment method, we use numerical quadrature to approximate the 


moments pk and the vectors Sk = u H pk, respectively defined in (14) and (19). 


3.4.1. Numerical approximation: trapezoid rule for moments 


Consider the equation (14) and assume that T has a 27r-periodic smooth 
parametrization: 


99 G C 1 (M, C), p{t + 27t) = p{t) Vf G M. 


Then, for k — 0,..., 2 m — 1, we have: 

hk = 2 “ f z k f(z)dz = ^ V(f) k f(v{t))v'{t)dt. 

Taking equidistant nodes tj — j — 0,..., N — 1, and using the trapezoid 
rule, we obtain the approximation: 

1 7V “ 1 
3 =0 

3.4.2. Numerical refinement: incorporating line search into Newton’s method 
Once an invariant pair has been numerically approximated, it can be 
refined using an iterative method such as Newton: this is, for instance, the 
strategy proposed in [ 3 ]. 

Newton’s method defines the correction (AX, AS) at each iteration as 

P(X, S) + DP (X , 5 ) (AA, AS) = 0 ( 21 ) 
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In this section, we show how to incorporate exact line searches into New¬ 
ton’s method for solving the invariant pair problem P(X, S ) = 0. Line 
searches are relatively inexpensive and improve the global convergence prop¬ 
erties of Newton’s method (see H). 


Algorithm 1. (Newton’s Method with Line Search) 

Input: initial approximation (Xq,So), tolerance e. 

Output: better approximation {X k ,S k ) to ([3]). 
step 1: Set k = 0 

step 2: If l|P( fe^ )llF < e- STOP 

step 3: Solve for {AX k ,ASk) the equation: 

HP (X ' S) (AX k , A S k ) = -P{X k , S k ) (22) 

step 4 : Find by exact line searches a t that minimizes the function: 

mm\\P(X+ tAX,S+ tAS)\\l (23) 

te[ 0,2] 

step 5: Update 

• X k+ \ = X k + t.AX k , S k+ 1 = Sk + tAS k . 

• k = k + 1 and go to step 2. 

Each iteration of a line search method computes a search direction d k and 
then decides how far to move along that direction. The iteration is given by 


■£fc+1 k T t k d k 


where the positive scalar t k is the step length. The success of a line search 
method depends on effective choices of both the direction d k and the step 
length t k (see 1321). A value t k = 1 gives the original Newton iteration. 

In our specific problem (|3]) , the direction d k is given by the solution 
(AX k ,AS k ) of the correction equation (22). The step length t k on each 
iteration is given by the solution of the minimization problem: 

pit) = \\P(X + tAX, S + f AS 1 ) |||. 


To facilitate the calculation, we consider the following equivalent representa¬ 
tion for ((3|) . A pair (X,S) G C nxfc x C kxk is an invariant pair if and only if 
satisfies the relation (see 0 ): 


1 

2 m 


l 


P(X)X{XI — S) dX 


0, 


(24) 
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where T C C is a closed contour with the spectrum of S in its interior. Using 
the formula for the total derivative of P at (. X , S ) in direction (AX, AS): 


DP (X|S) (AJf, AS) = J-<f P{ A) (AX + X(XI- S)- 1 AS) (AI-5) -1 dA, 

2m J r 

we have, at second order in ||AX|| and ||/A*S'||: 

P(X + tAX, S + tAS) = — I p(X)(X + tAX)(XI - S - tAS^dX = P(X, S) + m>P lx S) (AX, AS)+ 
2m Jr 1 ; 

— (f P( A) [AX + X(XI - S) _1 AS] (XI - S)- 1 AS(XI - S)~ 1 dX + 

_ 2m Jr 

— I P(X)AX(XI - S)- 1 AS(XI - S)~ 1 AS(XI - S) _1 dA 
2m Jr 


+ r 


+ r 


Recalling that Newton’s method defines (AX, AS 1 ) by (21), we have: 

P(X + tAX, s + tAS ) = (1 - t)P(X, S)+ 


+ r 


+ P 


— (f P(X) [AX + X(XI - 5)" 1 A5l (A/ - S)- 1 AS(A7 - 5) _1 dA 
2m 7r 

— I P(X)AX(XI - S)- 1 AS(XI - S)~ 1 AS(XI - Sp 1 dX 
2m Jr 


Thus, we have: 

p(t) =(1 - t) 2 \\P(X, 5)||| + P\\Af F + t 6 ||P|||. + 1 2 ( 1 - t)trace(P(X, S)*A + A*P(X , 5)) + 
+ 1 3 ( 1 - i)trace(P(X, 5)*P + P*P(X, 5)) + t 5 trace(XP + P*A) 

=(1 — t) 2 a + t 4 0 + t 6 <p + t 2 ( 1 — t)/3 + t 3 (l — t) 7 + t 5 r] 

=t 6 ip + t 5 r] + t 4 (0 - 7 ) + t 3 (7 — /?) + t 2 (a + /3) — 2crf + a (25) 

where: 

A= — l P( A) [AX + X(A/ - 5)" 1 A5l (A/ - S )" 1 A5(AJ - 5)" 1 dA, 

2m J r 

B = 2- / P(\)AX(XI - S)~ 1 AS(\I - S)~ l AS(XI - SP'dX, 

2m J r 

a = || P(X, 5)111., 9 = \\A\\ 2 f , <p = \\B\\ 2 f , r, = trac e(A*B + B* A), 
p = trace(P(X, S)*A + A*P(X, S)), 7 = trace(P(X, S)*B + B*P(X , 5)). 


Therefore, in each iteration of Algorithm [TJ solving (23) is equivalent to 
finding the minimum of the polynomial p(t) for t G [0, 2]. 


3.5. Numerical results 

In this section we compare two methods to refine approximate invariant 
pairs (X, S') G C nxk x C kxk : Newton’s method (N.M.) presented in [3] and 
Newton’s method with line search (N.M.L.S.), explained in Section 3.4.2[ 















We have implemented both methods in MATLAB and applied them to 
several problems taken from the NLEVP collection (see 0)- For each prob¬ 
lem, an initial invariant pair (X 0 , So) has first been approximated using the 
(block) moment method of Section 3.2 and approximating the moments /i, in 
(14) via the trapezoid rule discussed in Section 3.4.1 with N = 20 integration 
nodes. Moreover, F is chosen for each problem as the contour enclosing the 
k eigenvalues with largest condition number (computed using the MATLAB 
function polyeig). 

Table [ljshows that line search is generally effective in reducing the number 
of iterations and the overall computation time. 


N.M. N.M.L.S. 


Problem 

Deg P(A) 

n x k 

Ite 

Time 

Ite 

Time 

bicycle 

2 

2 x2 

23 

0.082 

16 

0.112 

butterfly 

4 

64 x 5 

67 

3.719 

22 

1.567 

cd player 

2 

60 x 6 

500 

N.C. 

19 

1.021 

closeddoop 

2 

2 x2 

8 

0.016 

7 

0.02 

damped beam 

2 

200 x 6 

28 

6.109 

4 

0.6037 

dirac 

2 

80 x 6 

500 

N.C. 

41 

1.965 

hospital 

2 

24 x 24 

53 

5.65 

51 

6.21 

metaLstrip 

2 

9x9 

500 

N.C. 

28 

0.589 

mobile-manipulator 

2 

5x2 

8 

0.014 

7 

0.030 

pdde stability 

2 

225 x 6 

29 

9.644 

16 

5.622 

planar .waveguide 

4 

129 x 6 

72 

11.148 

19 

3.682 

plasma drift 

3 

128 x 6 

69 

13.059 

26 

5.596 

power-plant 

2 

8 x8 

15 

0.34 

13 

0.39 

railtrack 

2 

1005 x 3 

32 

199.365 

28 

209.471 


Table 1: Comparison of results for classical Newton and Newton with line search. 

Figure [l] shows the convergence of the Newton’s method with line search, 
for the Dirac problem presented in Table [TJ Here we use as contour the circle 
of center C = —0.1 and radius R = 1.14, which contains the 6 eigenvalues 
with largest condition number. 


4. Matrix solvents 

In this section we study the matrix solvent problem as a particular case 
of the invariant pair problem, and we apply to solvents some results that we 
have obtained for invariant pairs. 
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Figure 1: Convergence of Dirac problem using Newton’s method with line search.This is 
a log-10 plot of the relative residual versus the number of iterations. 

Definition 6. A matrix S E C nxn is called a solvent for P(S) if satisfies the 
relation: 

P{S ) := A e S e + • • • + A 2 S 2 + AiS + 4) = 0. (26) 

A special case is, for £ = 2, the quadratic matrix equation Q(S) := 
A 2 S 2 + A\S + A 0 = 0, which has received considerable attention in the 
literature. For instance, in ra and EH the authors find formulations for 
the condition number and the backward error. They also propose functional 
iteration approaches based on Bernoulli’s method and Newton’s method with 
line search to compute the solution numerically. 

The relation between eigenvalues of P( A) and solvents is highlighted in 
[25]: a corollary of the generalized Bezout theorem states that if 5 is a solvent 
of P(S) then: 

P(A) = L(A)(A/-S), 

where L( A) is a matrix polynomial of degree £ — 1. Then any eigenpair of 
the solvent S is an eigenpair of P{ A). 

4-1. Condition number and backward error 

An analysis and a computable formulation for the condition number and 
backward error of the quadratic matrix equation Q(S) = 0 can be found in 
pH! and [21]. 
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Here we give explicit expressions for the condition number and backward 
error for the general matrix solvent problem P(S) = 0. We follow the ideas 
presented in [20] . [2Tj and [36] . 


4-1.1. Condition number 

We perform here a similar analysis as we did in Section 2.1.1 
A normwise condition number of the solvent S can be defined by: 

k(S) = limsup {- ^ : (P + A P)(S + AS) = 0, \\AAi\\ F < eai,i = 0 : (27) 

e->0 l 6 I P 11 _F 


e 

where AP(A) = A*A Ai. The on are nonnegative weights; in particular, 

i =o 

AAi can be forced to zero by setting a* = 0. 


Theorem 8. The normwise condition number of the solvent S is given by: 

B~ s 1 B a 

<s) = 1ST' 

where 

i j -1 

= [a e (S e ) T ®In a t -i (S^Y ® In ■■■ a 0 /„*] 

3=1 i=0 

The proof of Theorem [8] follows from the proof of Theorem [2j by taking 
A A" = 0, X = / and noting that the matrix S has size n x n in the matrix 
solvent problem. 


4-1.2. Backward error 

Let ci!j, for i = 0, be nonnegative weights as in Section 4.1.1 The 
backward error of an approximate solution T to (26) can be defined as: 

V(T) = min{e : (P + AP)(T) = 0, ||AAj|| F < ea^i = 0,... ,£} (28) 

We proceed as in Section |2.1.2| and we obtain bounds for the backward error 
of P(S): 

j/(T) > 1'^™- 


ri(T) < 


( a I W^Wf s-l «?II t IIf + «o) 1/2 ’ 


(aja min (T £ ) 2 + ''' + «?a min (T)2 + a g) 1/2 ’ 

where a m ; n denotes the smallest singular value, that by assumption is nonzero. 
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5. Computation of solvents 

Motivated by applications to differential equations [10] . we study an ap¬ 
proach to the computation of solvents based on the moment method, by 
specializing the results presented in Section |3j 

Let us recall some results that will be needed later. The next result is a 
generalization of a theorem presented in EH which gives information about 
the number of solvents of P(S). 

Theorem 9. Suppose P( A) has p distinct eigenvalues {X i} p i=1 , with n < p < 
in, and that the corresponding set of p eigenvectors {vi}? =1 satisfies the Haar 
condition (every subset of n of them is linearly independent). Then there are 

at least different solvents of P(\), and exactly this many if p — in, 

which are given by 

S = Wdiag(m)W-\ W — [wi ■■■ w n ] , 

where the eigenpairs (pi, uy)” =1 are chosen from among the eigenpairs (A*, u«)f =1 
ofP. 

Note that if we have that p — n in Theorem [9j the distinctness of the 
eigenvalues is not needed, and then we have a sufficient condition for the 
existence of a solvent. 


Corollary 2. If P( A) has n linearly independent eigenvectors V\,V 2 ,... ,v n 
then P(S) has a solvent. 


An example which illustrates this last 
the quadratic matrix solvent problem (see 


result is the following. 

ra. ed) 


Consider 


Q(S) = S 2 + 




12 

14 


Q( A) has eigenpairs: ^1, ^ 
example, the complete set o: 


2 , 


solvents is: 



and ( 4, 


1 

1 


. For this 


1 0 


1 

2 


3 

0 


1 

3 

and 

0 2 

J 

0 

3 

1 

1 

2 

1 

0 

4 
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Note that we cannot construct a solvent whose eigenvalues are 3 and 4 be¬ 
cause the associated eigenvectors are linearly dependent. 

Our approach to compute matrix solvents is based on the relation between 
the matrix solvent problem (26) and the invariant pair problem ([3]). We state 
this in the following result. 

Theorem 10. Let P( A) be a n x n matrix polynomial and consider an in¬ 
variant pair ( Y,T ) G C nxfc x C kxk of P( A). If the matrix Y has size n x n, 
i.e. k = n, and is invertible, then S = YTY 1 satisfies equation (26), i.e., 
S is a matrix solvent of P( A). 

Proof. As ( Y,T) G C nxn x C nxn is an invariant pair of P(A), we have: 

A e YT l + ■ • • + A 2 YT 2 + AiYT + A 0 Y = 0. 

Since Y is invertible, we can post-multiply by Y~ x . Then we get: 

A(YT l Y~ x + • • • + A 2 YT 2 Y ~ 1 + AiYTY -1 + A 0 = 0 
A,S l + • • • + A 2 P 2 + AtS + A 0 = 0. 

Therefore, S' is a matrix solvent of P( A). □ 

5.1. Numerical refinement of solvents 

As pointed out before, the use of Newton’s method incorporating line 
searches to find solvents is not new. For instance, in [2T], [27] this approach 
is used to approximate solvents for the quadratic matrix equation. Here we 
apply this method to approximate solvents for the general matrix solvent 


problem P(S) = 0 and follow the ideas of Section 3.4.2 


The application of Newton’s method with line search to find solvents is 
based on the following steps for the fc-th iteration. 

• Solve for ASk the equation: H)Ps(ASk) = — P(Sk ) 

• Find by exact line searches a tk that minimizes the function: 

min || P(S k + t k AS k )\\ 2 F 

tk 6 [ 0 , 2 ] 


• Update S k+ i = S k + t k AS k . 
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The step length on each iteration is given by the solution of the minimiza¬ 
tion problem: 

p{t k ) = \\P{S k + t k AS k )\\ 2 F , 


As we did in Section 3.4.2| for the invariant pair problem, we use an equivalent 
contour integral representation for (26). A matrix S £ C nxn j s a solvent if 
and only if satisfies the relation: 


2- (f P(A)(A/ - S)~ 1 d\ = 0, 

2m J p 

for any closed contour T C C with the spectrum of S in its interior. 
Using the formula for the total derivative of P at S in direction AS: 


BP s {AS) 


1 

2m 


P(A)(M - 5)-'AS(A/ - 


(29) 


we obtain: 


P(S + tAS) 


p(S) + mp s (AS) + t 2 


i 

2iri 


i 


P(A)(A/ - S )- 1 AS(\I - S)~ 1 AS(XI - Sy 1 d\ 


Recalling that Newton’s method defines AS by 


P(S) + BP s (AS) = 0 


then we have: 


P(S + tAS ) 


(1 -t)P(S) + t 2 


— (f p(x)(xi - s)~ 1 as(xi - sy'Astxi - sy'dx 

2m J T 


Thus, we obtain: 


p{t) = (1 - t) 2 \\P(S)\\ 2 F + t 4 \\A\\ 2 F + t 2 ( 1 - t)tiace(P(S)*A + A*P(S)) 
= (1 — t) 2 a + t A 0 + t 2 (l — t)/3 = 

= t 4 9 — t 3 /3 + t 2 (a + j3) — 2 at + a 


where: 

A — — / P(A)(AJ - S)~ 1 AS(XI - S)~ 1 AS(XI - 5) _1 dA, 

27u 

a= ||P(S , )|£, e=\\A\\ 2 F , (3 = trace(P(S')*A + A*P(S)). 

Therefore, in each iteration one finds the minimum of the polynomial p(t) 
for t £ [0, 2], 
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6. Solvents and triangularized matrix polynomials 

Motivated by the results in [35] and [38], where the authors analyze a 
method for triangnlarizing the matrix polynomial P(A), we aim here to study 
the relation between solvents of general and of triangularized matrix polyno¬ 
mials. 


6.1. Triangnlarizing matrix polynomials 

For any algebraically closed field F, any matrix polynomial P(X) E F[A] nxm , 
with n < m, can be reduced to triangular form via unimodular transforma¬ 
tions, preserving the degree and the finite and infinite elementary divisors 
[35], [38]. 

Theorem 11. [35 j/ For an algebraically closed field F, any P(A) E F[A] nxm 
with n < m is triangularizable. 

What can we say about solvents for a given matrix polynomial P(A) and 
for the associated triangularized polynomial? A partial answer will be given 
in Theorem [13] 


Theorem 12. For any faxfn monic linearization XI —A of P( A) e C[A] nxn 
with nonsingular leading coefficient, there exists U E C nxfe with orthogonal 


columns such that M 


U 

UA 


is nonsingular and XI — MAM 1 is a 


UA *- 1 

linearization for the polynomial T (A) = X e I + A* -1 T£_i 4 -bA 2 X 2 + ATi-|-To, 

which is upper triangular and equivalent to P(A). 


Theorem 12 is a straightforward generalization of a result found in 
Note that, for the time being, we assume that the leading coefficient An is 
nonsingular. 


Theorem 13. Let P(A) be a n x n matrix polynomial and consider the 
linearization: 

'0 /„ 0 0 

0 0 I n ■■■ 0 


A = 


0 0 0 ■ ■ ■ I n 

—A 0 —A\ — A 2 ■ ■ ■ —An -1 
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M 


-i 


St 1 


Let M be as in Theorem 12 and let Yj be the first n x n block of the matrix 

In 

St 


If Y\ is nonsingular and S t is a solvent for the triangularized 


problem, i.e., T(S t ) = 0, then S = Y x S f Y x 1 is a solvent for P(S). 
Proof. Note that: 


'O' 


St - St 


0 

_ 

s 2 t - 4 2 


0 


-To -T x S t - T 2 Sf - • 

• - X 


0 I n 0 

0 0 I n 

0 0 0 

-T 0 —T-i -T 2 


0 

0 

In 

-Tt-i 



' In ' 


[41 


St 


4 2 


St 1 . 


si 



0 

4 

0 

0 


‘ T " 


f 


0 

0 

In - 

0 


■*■71 

s t 


■*■71 

s t 


0 

0 

0 

In 

MM- 1 

1 

.,. 7 

^ -+o 
c °. 

1 

nt- 1 



-Ti 

-T 2 ■■ 

• -4-i. 



L4 J 


St = 


(Hi) 


0 

In 

0 

0 


T 

r 

0 

0 

4 •• 

0 


-*n 

St 


O • • • 

0 

0 

4 

M -1 

i 

Co 

1 ' ■ ■ - 

1_ 

-M~ x 

i- 

i 

o 

-Ai 

— a 2 ■ • 

i 

7 

1 


- 


In 

St 


S 


Since M 


-l 


In 

St 


S 


r£-l 


has size In x n , let us partition it as 


C nxn for i = 1,..., I. Then: 

' o 4 o 
o o /„ 

ooo 

—Aq —A\ —A 2 


<t-i 

X 

4 

Y 


S t . 


, where 4 € 


0 

0 


— An-\ 



v 


vr 


4 

_ 

y 2 


Y e _ 


_4 


S t = 


36 







































Then we have: 


Y i = Y 1 S^\ for i — 2,... ,£‘ 

-AqYi A i Y‘2 - At^Yi - Y e S t = 0. 


Substituting equations (30) in (31) we obtain: 


0 — Y\S t + At_\Y\S t + • • • + A{Y\St + AqY\. 

If Y\ is invertible we have: 

0 = Y 1 SfY 1 ~ 1 + Ag-iYiSf^Yf 1 + • • • + A^StYf 1 + H 0 . 
Taking S = Y\S f Yfi l we have: 

0 = S l + A^S 1 - 1 + • • • + A 2 S 2 + AxS + Ao := P(S). 
Then S = YiS t Y^ 1 is a solvent for P(S). 


(30) 

(31) 


□ 


6.2. Example: A problem with an infinite number of solvents 

What happens to the ideas outlined above when working on problems 
with an infinite number of solvents? Here is an example taken from [30| . 
Consider the matrix polynomial: 


P{ A) = A 2 / + A 


r-7 

-2 

-2' 


r 13 

9 

7 1 

3 

-203 

8 

+ 

-21 

294 

-36 

31 

-13 

31 

-40 

31 

-231 

31 

60 

31 

183 

31 

435 

L 31 

31 

31 J 


L 31 

31 

31 J 


Triangularizing T(A) we find: 



"(A — 3)(A — 

4) (A- 

3) 

0 


T(A) = 

0 

(A- 

3) 2 

1 

= 


0 

0 


(A — 4) 2 _ 




"-7 

1 0 ' 


'12 -3 

O' 

= A 2 / + A 

0 

-6 0 

+ 

0 9 

1 



0 

0 -8 


0 0 

16 


A 2 / 2 + AT X + T 0 = 


Now, suppose that the solvent S t £ C 3x3 of the triangularized problem is in 
upper triangular form, i.e.: 


S t = 


Xu 

0 

0 


Xl2 

X22 

0 


£13 
£23 , 

£33 
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then we have: 

T(S t ) =S t 2 +T 1 S t +T 0 = 

"(^ll — 3)(xn — 4) X 2 2 — 7^12 + xnx 12 + x 12 x 2 2 ~ 3 x 2 3 — 7xi 3 + xnaj 13 + xi 2 x 2 3 + X 13 X 33 ' 

= 0 (x 2 2 — 3 ) 2 X22X23 ~ 6x 2 3 + ^23^33 + 1 

0 0 (x 33 - 4 ) 2 

In the task of solving the problem T(S t ) = 0, we see that: Xn = 3 or £ n = 4, 
X 22 = 3 and x 33 = 4. Then we have two cases: 

I. If X\\ = 3, X 22 = 3 and £33 = 4: 

Then we find that £23 = — 1, £12 = 0 and £12 = — 1, which is a contra¬ 
diction. In this case there is no solution and then we can’t construct a 
solvent. 

II. If £n = 4, £ 2 2 = 3 and £ 33 = 4: 

Then we find that £23 = — 1 and £13 = £12 + 1. In this case the solvent 
St has the form: 


4 £12 £12 + 1 


O 

_1 


"0 1 1' 

0 3 -1 

= 

0 3-1 

+ X 12 

0 0 0 

O 

O 


O 

O 


0 0 0 


for £ 12 € C. 

Thus T(A) has an infinite number of solvents and the same holds for P( A). 

7. Conclusion 

I 11 this paper we have explored several questions related to invariant pairs 
and solvents of matrix polynomials. In particular, preliminary results on the 
use of scalar or block moment Hankel pencils to compute invariant pairs and 
solvents suggest that this approach may present several points of interest. A 
more detailed analysis, along with the design and development of effective 
algorithms and extensive numerical tests, will be the topic of future work. 
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